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Abstract 

We simulate Grover's algorithm in a classical computer by means 
of a stochastic method using the Hubbard-Stratonovich decomposition 
of n-qubit gates into one-qubit gates integrated over auxiliary fields. 
The problem reduces to finding the fixed points of the associated sys- 
tem of Langevin differential equations. The equations are obtained 
automatically for any number of qubits by employing a computer al- 
gebra program. We present the numerical results of the simulation for 
a small search space. 

1 Introduction 

Grover's and Shor's algorithms ^1 12] are landmarks in the development of 
Quantum Computation since the foundation of this area in the 1980's [HlllllH]- 
Shor's algorithm provides an exponential speed up over currently known clas- 
sical factorization algorithms. Although a classical polynomial factorization 
algorithm might be found, it is very unlikely that its complexity would be 
as low as Shor's 0{n^). Compare with the polynomial 0{n^'^ poly(logn)) al- 
gorithm for the much easier problem of primality test recently developed by 
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Agrawal, Kayal, and Saxena [U]. Therefore, it is unlikely that a polynomial 
factorization algorithm, if it exists, would have complexity under 0(n^^). 

Shor's algorithm is the best candidate to be simulated in classical com- 
puters, since the price paid by the simulation may still be low compared to 
the exponential timing of classical algorithms. Unfortunately, the techniques 
used in this paper do not seem to be suitable for Shor's algorithm [7j. 

Grover's algorithm plays an important role in Quantum Computation, 
since it provides a proof that quantum computers are faster than classical 
ones for unstructured database searching. It has complexity 0{^/N) while 
the best classical algorithm has complexity 0{N), where is the database 
size. Grover's algorithm is optimal using oracle searching jHl [HI, therefore 
the quadratic speed up is the best one can achieve in this case. The hope 
in simulating it in classical computers with a better efficiency compared to 
the classical algorithm is very low. Despite this fact, Grover's algorithm 
provides a laboratory to test simulation techniques that might shed new 
light in distinguishing the way classical and quantum computers work. 

In a seminal paper, Feynman |3] argues that classical computers can sim- 
ulate quantum ones only with an exponential slow down. It is easy to see 
that an exact simulation uses either an exponential amount of memory or 
takes an exponential time to accomplish generic tasks that could be done 
by a quantum computer. A generic quantum state of n qubits has 2" com- 
plex amplitudes which must be taken into account by the classical computer. 
An alternative option is a stochastic simulation by some "non-deterministic 
computer". The central issue is the local nature of classical computers. A 
classical computer of the same physical size of a quantum computer work- 
ing under local classical physical laws does obey the Bell inequality jTU] and 
cannot reproduce the same results of a computer which follows the quantum 
mechanical laws. 

Feynman's arguments, on the other hand, are not quantitative. The ques- 
tion is: what is exactly the loss introduced by the stochastic simulation, in 
such a way that the result given by the quantum computer is still reproduced? 

An interesting stochastic simulation technique in the quantum computa- 
tion context was introduced by Cerf and Koonin [7j. The quantum computer 
is viewed as a many-body dynamics which can be reduced to the time evolu- 
tion of single qubits integrated in auxiliary fields. They used the Hubbard- 
Stratonovich [TTl [T2] representation to obtain an expression for two-qubit 
gates in terms of two one-qubit gates with two auxiliary fields . They used a 
Monte Carlo Method fact the Langevin equations with equivalent 
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solutions [ini Cni) to reduce the exponential size of the associated Hilbert 
space to an equivalent coupled system of differential equations with a poly- 
nomial number of auxiliary fields. 

It is known that a general quantum circuit can be decomposed in terms 
of one and two-qubit gates, called universal gates [T71 HHj. On the other 
hand, it is also known that this decomposition may use an exponentially 
large number of universal gates ^H] (and auxiliary fields). In that case, the 
gain provided by the stochastic method would be lost in the decomposition 
in terms of the universal gates. 

In this paper, we generalize Cerf and Koonin's method to general n-qubit 
gates, avoiding the decomposition of n-qubit gates into universal gates. The 
generalization is straightforward and welcome in simulating Grover's algo- 
rithm for a general number of qubits. We apply the method for Grover's 
algorithm, and make use of a Maple program that generates automatically 
the associated system of Langevin equations. The equations are solved nu- 
merically (in Maple and in C) giving the relaxation values of the auxiliary 
fields which correspond to the fixed points of the Langevin dynamics. We 
present the numerical results for the simplest case of Grover's algorithm and 
discuss some problems we are facing with the method. 



2 Decomposition in terms of one-qubit gates 

Consider a quantum circuit of g gates (f/^, A; = 1, g) in a n-qubit quantum 
computer. Each gate Uk acts on two or more qubits {jk^\3k'\ it'^'')- The 
qubit indices within a gate Uk are put in parentheses. They do not coincide 
necessarily with the overall qubit indices, which will be denoted in square 
brackets. 

In order to clarify the notation consider the following example: if Uk is 
a two-qubit gate acting on the third and fifth qubits, then 77.^ = 2, = 3 
and j^j^^ = 5. The index s of j'l^!'^ runs from 1 to n^. The actual values of 
j^^^ are in increasing order but not necessarily consecutive with respect to 
the overall qubits, as shown in the example. One-qubit gates will be trivially 
introduced in the next section, so they will not be considered for a while. 

The computation as a whole is performed by the unitary operator U given 
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by 

1 

U = l[Uk = Ug ...U, (1) 

k=g 

where the product is in reverse order, so the operators act in ascending 
order from left to right on input kets. Uk is a general n^-qubit gate, with 
2 < Tifc < n. We assume that Uk can be decomposed as an exponential of a 
tensor product of one-qubit gates, in the form 

where A^*^ acts only on qubit This equation generalizes eq.(2) of Cerf 
and Koonin's paper [7^. 

The Hubbard-Stratonovich decomposition fn\ IT^ for a Ti^-qubit gate 
follows from 

g-*a,AW® ... ^a(=' _ r ^^(s) g-ia^AW® ... ^A^^-^)®^^'/'^' ^ 

J — oo 

where I^'^ is the identity matrix acting on the qubit and from the fol- 
lowing representation for the Dirac delta function 

- r dri'^ a, ,-...r<»'^'0...<-'^(Al='-4=)/<^)). (4) 

Applying decomposition (jS)) and (@)) recursively n^, — 1 times in eq.Q, we 
introduce real auxiliary fields a^^^ and r^'*'', 2 < s < n^. For simplicity, we 
will drop the identity matrices from now on. The result is the following 
expression: 

Uk = frfaf drf ...cia^)rfr^) ^ 
Substituting ^7^ into eq.((H) yields 



where 



vf'fe) = k!' Vf.^^. 4"*') -e-""' 



I'fV") = e-">'"l'><', 2<m<nt. (7) 

CTfc denotes all the alf^ for a given fc, and the superscript (s) indicates in which 
qubit A^f^^ acts. 

Up to this point, we have been considering U as a. composition of g gates. 
Using the Hubbard- Stratonovich representation, each gate was decomposed 
into one-qubit gates, paying the price of introducing 2(nfc — 1) scalar fields 
and respective integrals. Now we focus on the world line of each qubit and 
multiply (horizontally) all one-qubit gates acting on that qubit. We end up 
with n one-qubit gates which are tensored out and integrated on the auxiliary 
fields, yielding the following expression for U 

/n 
Da Dt ^"^^ '^l^'^^' (g) /J^ (a, r) (8) 

1=1 

where square brackets denote an overall qubit label (as opposed to internal 
gate labels jl^^), with 

U^H^,r) = llu!\a,T), (9) 

k=g 

(here a and r denote all of them) and 

y«(a,) if/ = ji 



(10) 



otherwise. 



In the remaining of this section, we follow closely Cerf and Koonin's 
approach [7]. The last step of the computation is the measurement of the 
first m qubits, m < n. We assume that the observable is the direct product 
of m one-qubit observables. 



C) = (g)CW. (11) 



1=1 



/^From the remaining qubits, p < n — m have prescribed values tti, 7r2, . . . , vr^, 
so we define the projector V as 



p 



r = l[Vk, (12) 

k=l 

where Vk = Ivr^) {TTk\. The expectation value of the observable O is given by 

- '"-■■""if::?^^!"-^"?' (13) 



Using eqs.dHI), (US) and (HHI) we get 

^ / Da Dt Da' Dr' exp(-25(a, r, a^ r')) Oja, r, a', r') 
^ ' J Da Dt Da' Dt' exp{-iS{a,T,a',T')) ^ ' 

where 

S'(cr,T,cr',T ) = 

fc=l \s=2 
n 

+1 J2 In (Oi I U^^'^ [a', r')Pt'l t/W (a, r) |0,) (15) 
(=1 

and 

0{a, T, a', r) = I I ^ , Z r.mrm , W 



C — 

k=l \s=2 



(0/1 


mH(a',T')CWpWt/W(^, 






\ Wli] {a', r')PmV] {a, t) 


m 



where O^^^ = I for / > m. 

The most promising method to calculate (O) is by using the associated 
Langevin equations given by 

(s) 

where t is a simulation time, is a real Gaussian white noise and the scalar 
fields have been extended to the complex plane. The stochastic estimate of 
(O) is 

{0)^^J^'^^dtOiait),a'it)) (18) 
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where a and a' are solutions of eq. lfTTj) . and T represents a sufficiently large 
simulation time. In the next section, we will obtain that a and a' converge 
to some final limiting value, so we will use the formula 

(0)^0(afinai,4,,i). (19) 

3 Simulation of Grover's algorithm 

Grover's algorithm allows us to search an element in an unstructured 
database with elements (suppose that N = 2"') using C(a//V) steps. The 
best classical algorithm uses 0{N) steps. This quadratic speed-up is one of 
the greatest successes of quantum computation so far. 

Grover's algorithm has two registers, the first one with n qubits and the 
second one with one qubit. It starts by preparing a superposition of all 
computational basis states with same amplitude. An oracle is used to probe 
the database (a quantum memory is assumed to exist in some form), and 
it changes the sign of the amplitude of the state which corresponds to the 
numerical index of the searched element. In order to simulate the oracle 
action, we use a t?, + 1 qubits generalized Toffoli gate, which marks the state 

— 1) (we are assuming that — 1 is the searched element). 

The state |A^ — 1) has n ones in binary notation, and that is why the 
oracle just changes the sign of the amplitude of that state. Any other state 
of the computational basis could be marked as well, introducing two X gates 
in the oracle, for each zero in the binary representation of the state's label. 
The X gates are placed symmetrically with respect to the generalized Toffoli 
gate, in the qubits corresponding to the zeroes in the binary representation. 

Fig. ^ shows the Grover operator in terms of two generalized Toffoli gates 
and one-qubit gates. This decomposition is enough to apply the method of 
the previous section. Note that the method used in [7] would require the 
decomposition of the generalized Toffoli gate into two- and one-qubit gates. 
This would introduce more scalar fields (a's and r's) and lead to cumbersome 
calculations. For example, for A^ = 4, the method of Cerf and Koonin would 
use 28 scalar fields, while ours uses 12. 

A non-trivial part of our procedure is the calculation of the decompo- 
sition (0) for a general n-qubit gate. Fortunately, it is easy to obtain the 
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Figure 1: Grover operator. 
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decomposition of the generalized Toffoli gate: 



where 



/ 1 



a 



1 

1 oy 



-iayl(i)®...®A("+i) 



2' 



1 - cr.r 



2 

1 -1 
-1 1 




1 



(20) 



f211 



The next step is the calculation of {/['^(cr, r) given by eq.© for I = 
1, . . . ,nfc. Recall that [/[^'(cr, r) is the ordered composition of all one-qubit 
gates for the qubit including the two gates that come from the decompo- 
sition (j^Uj) . one for each generalized Toffoli gate of Fig. |21 
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Figure 2: Full Grover algorithm. 



Since the Grover operator G must be apphed ko times, we have 

{UUoUh^^^r)) H for/ = l,2,...,n 

XH(uUoUt^'\^,r))HX foTl = n+l. 



(22) 



Examining fig^ we see that there are four distinct kinds of [/t'^(cr, r), which 
we address now. For / = 1, [/H(a,T) must be treated separately, since from 
eq. (fTO|l we see that for / = ji (here ji = 1) V^^^ depends only on the a fields 



while V^^^ depends on tI'^\ So 



up = HX 



e" 



/2) (n) 



XH 



e 



(n + 1) 



(23) 



For I 



, n — 1 the gate configuration is similar, then 
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.(0 



(24) 



For I = n, Grover operator has two extra Hadamard and a X gate. Then 



= -HXH 



Finally, for I = n + 1, 




l-e 
1 + e 



-2ia T, 



(n) 



-lia T,- 



(n) 



HXH 



1 
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An) 



n+l] 



T (n + l) 

1 _ e~^*"^2fe-i 



1 -2ja r<"+'' 

1 + 2fe-i 



Next we can insert the latter results in eq.(j22I) and obtain t/W. We note 
that in eq. (fTB|) . |0;) is represented as ^ . So we get for example, for 
observable Oqo = (|0i) (Oi|) ® (IO2) (O2I), 



Coo (a, r, (r',r') 



n 



t/t'"i,i(a',r')f/[']i,i(a,T) 



n 



m^,i(a',rO^Wi,i(a,r) + f/tWi,2K,r')^W2,i(^,r)' 



(27) 



1=1 



Similarly, in each observable we choose to test, we need to calculate only the 
relevant elements of f/t'l 

In the Appendix, we describe a Maple^ program that calculates the phys- 
ical quantities of this section starting from this point. Given the program 
uses eqs. (I2SI) to (|2ni) to calculate [/W(a,T) and [/^ (cr',r'), and then obtains 
0{a^ T, a', t') for a given observable constructed similarly to the above. Us- 
ing eqs.(|15p and (jl7|) . the program obtains the Langevin equations for the 
simulation of Grover's algorithm with an arbitrary number of qubits. 

The simplest case is a search space of = 4 elements, for which the auxil- 
iary fields are (t^\ o"(^\ '^2'\ and the respective primed versions. 
We give below the first Langevin equation for this case. The other equations 
of the system of partial differential equations have a format similar to this 
one: 



This system of equations is discretized and solved numerically. After trying 
different values for the discretization parameter and initial conditions, we plot 
each scalar field a\^^ as a function of t in order to determine the convergence 



da 



dt 




(28) 



^Maple Waterloo Software, Inc. See http://' 



www.niaplesoft.com 
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value, to be substituted in the previously obtained formula (j27j), for the 
observable (Cqo) and in the equivalent formulas for (Oqi), (Cio) and {On), 
corresponding to the observables 

Oi, = \ij){ij\ (29) 

for i,j = 0, 1. Since we have chosen the oracle that changes the sign of the 
state 1 11) |— ), we expect to obtain (On) close to 1, and the remaining ones 
close to 0. 

Unfortunately, we are faced with technical difficulties. First, there are 
some fields that do not converge to a definite value, at least to the extent 
of our simulation. In our case, this issue was circumvented by imposing the 
constraint that total probability equals one, but in the general case the use 
of such an artifact would not be desired. The second problem are fields 
that seem to have a logarithmic behaviour. We have tried to solve this by 
fitting the logarithmic behaviour by a log function, then inserting it in the 
expressions and taking the limit as t ^ oo. Even then, contrary to our 
expectation, we get the following results 

(Coo) = .28 

{Ooi) = .24 (30) 
(0io) = .24 
(On) = .21 

It is not clear to us whether we are making some mistake or this result shows 
that the method does not work with Grover's algorithm. Any comment on 
this issue is very welcome. 
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Appendix 



This appendix describes a Maple program for release 6 or higher that cal- 
culates the Langevin equations for the simulation of Grover's algorithm, as 
described in Section IHl The number of qubits is fixed in the beginning of the 
session. The lines beginning with the prompt > are Maple commands. We 
make some comments about each group of commands. 

1. Starting a new session, defining matrices H, X and Z, and setting the 
number of qubits. 

> restart; 

> with (Linear Algebra) : 

> H := l/sqrt(2)*Matrix([[l,l] , [1,-1]]) : 

> X := Matrix([[0,l] , [1,0]]) : 

> Z := MatrixC [[1,0] , [0,-1]]) : 

> n := 2: 

2. Function kO finds the number of times the Grover operator is applied. 

> kO := procO 

> local theta; 

> theta := 2*evalf (arccos (sqrt (l-l/(2~n) ) ) ) ; 

> eval(round(Pi/(2*theta) - 1/2)) 

> end proc: 

3. The next two procedures calculate [/['^(cr, r) given by eq.(j221) and uf\ak, Tk) 
given by eqs. (j23ll2'(3|) . They are not used directly, as we see ahead. 

> Ul := proc(l) 

> if Kn+1 then 

> ' . ' (seq(Ul_k(l,kO()-s) ,s=0. .kO()-l)) .H 

> elif l=n+l then 

> X.H. ' . ' (seq(Ul_k(n+l,kO()-s) ,s=0. .kO()-l)) .H.X 

> else error (' expecting l<=n+l, got %1',1) 

> end if 

> end proc: 
> 

> Ul_k := proc(l,k) 
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> if 1=1 then 

> H . X . Matrixl (muKsigma [s , 2*k] [t] , s=2 . . n) ) . X . H . 

> Matrixl (mul(sigma[s,2*k-l] [t] ,s=2. .n+1)) 

> elif Kl and Kn then 

> H . X . Matrixl (tau [1 , 2*k] [t] ) . X . H . Matrixl (tau [1 , 2*k-l] [t] ) 

> elif l=n then 

> Z.Matrix2(tau[n,2*k] [t] ) .Z. Matrixl (tau [n,2*k-l] [t] ) 

> elif l=n+l then 

> Matrix2(tau[n+l,2*k-l] [t] ) 

> end if 



> end proc: 

4. The next two procedures are auxiliary functions for calculating the ma- 
trices of eqs.dHESI)- 

> Matrixl := proc(x) 

> Matrix([[l,0] , [0,exp(-I*(Pi/2)*x)]] ) 

> end proc: 
> 

> Matrix2 := proc(x) 

> MatrixC [[l/2+l/2*exp(-I*Pi*x) , -l/2*exp(-I*Pi*x)+l/2] , 

> [-l/2*exp(-I*Pi*x)+l/2, l/2+l/2*exp(-I*Pi*x)]] ) 

> end proc: 

5. The next two functions calculate [/'''(a, r) and its hermitian conjugate, 
automatically applying simplifying functions. An user should use these func- 
tions since they return the simplified result. 

> U := 1 -> Map(simplify,Ul(l)) : 

> Udagger := 1 -> Map (simplify, subs(l=-l, sigma=sigmal , 

> tau=taul, Transpose(Ul(l)))) : 

6. The next command calculates the action S given by eq. ()15|) . taking |0) as 
the input state for all qubits. 

> S := - (Pi/2)*add( 

> add(sigma[s ,k] [t] *tau [s ,k] [t] - 

> sigmal [s,k] [t] *taul [s,k] [t] , s=2. .n+2*frac(k/2)) , 

> k=l. .kOO) + 

> I* add( ln( simplify ( Udagger (1) [1 , 1] *U(1) [1 , 1] + 

> Udagger(l) [1,2]*U(1) [2,1] ) ), 1=1.. n+1): 
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7. The next commands calculate the Langevin equation given in eg . (11711 . 

> dS := map (f actor, diff(S,sigma [2,1] [t] )) : 

> Diff (sigma[2,l] [t] ,t) = map(x->x*(- 1/2), dS) + eta(t) ; 
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